Solution of an associating lattice gas model with density anomaly on a Husimi lattice 
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We study a model of a lattice gas with orientational degrees of freedom which resemble the 
formation of hydrogen bonds between the molecules. In this model, which is the simplified version 
of the Henriques-Barbosa model, no distinction is made between donors and acceptors in the bonding 
arms. We solve the model in the grand-canonical ensemble on a Ifusimi lattice built with hexagonal 
plaquettes with a central site. The ground-state of the model, which was originally defined on 
the triangular lattice, is exactly reproduced by the solution on this Husimi lattice. In the phase 
diagram, one gas and two liquid (high density-HDL and low density-LDL) phases are present. All 
phase transitions (GAS-LDL, GAS-HDL, and LDL-HDL) are discontinuous, and the three phases 
coexist at a triple point. A line of temperatures of maximum density (TMD) in the isobars is found 
in the metastable GAS phase, as well as another line of temperatures of minimum density (TmD) 
appears in the LDL phase, part of it in the stable region and another in the metastable region of this 
phase. These findings are at variance with simulational results for the same model on the triangular 
lattice, which suggested a phase diagram with two critical points. However, our results show very 
good quantitative agreement with the simulations, both for the coexistence loci and the densities of 
particles and of hydrogen bonds. We discuss the comparison of the simulations with our results. 

PACS numbers: 05.50.-f q,61.20.Gy,65.20.-w 



I. INTRODUCTION 



The introduction of orientational degrees of freedom 
in lattice gas models may result in rich phase diagrams. 
As an example, we may mention the study of lattice gas 
models with direction dependent interactions which were 
found to exhibit closed loop coexistence curves such 
as the ones found in solutions of glycerol with guaia- 
col 0, TO-toluidine [31 , and ethylbenzylamine 0, which 
exhibit a nearly symmetric coexistence loop with both 
an upper and a lower critical solution temperature. It 
was suggested by Hirschfelder, Stevenson, and Eyring 
that the low-temperature critical point might be due 
to a highly directional short-range interaction, such as a 
hydrogen bond: while at low temperature the ordering 
of these interactions lowers the energy of solution, with 
rising temperature this ordering is decreased and phase 
separation occurs. These suggestion was followed in the 
model proposed by Barker and Fock some time later 
and the solution of this model in the quasi-chemical ap- 
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proximation actually displays a closed coexistence loop. 
A simplified version of the model defined on a conve- 
niently decorated simple cubic lattice, may be mapped on 
the three-dimensional Ising model and thus much precise 
information is known about its thermodynamic behavior 
_6!j . When the directionality of the part of the interactions 
in the model due to the hydrogen bonds is increased, the 
results are closer to the experimental data for the mix- 
tures cited above, although the correspondence to the 
Ising model is lost 

In water, the ordering of hydrogen bonds is supposed 
to be important in determining the unusual thermody- 
namic and dynamic behavior, including the possible ex- 
istence of an experimentally unaccessible liquid-liquid 
phase transition Q ■ Liquid- liquid phase transitions were 
originally found by Monte Carlo simulations of realistic 
liquid water models with atomic details [1, Q , but they 
were already observed experimentally in systems such as 
phosph orus [l0| . triphenyl phosphite [111] and n-butanol 
[l3,[l3|- Tetrahedral liquids, such as silica and water, also 
present thermodynamic and dynamic anomalies which 
can possibly be related to the second critical point (SCP) 
associated with these transitions [ll-[l3. Among these 
anomalous features, we note the increase of density with 
temperature that happens in liquid water at tempera- 
tures below 4°C and the apparent divergent behavior of 
thermodynamic response functions with decreasing tem- 
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peratures towards the deep super-cooled liquid, at atmo- 
spheric pressures Q- 

Several lattice models with orientational interactions, 
usually called network-forming fluids or associating lat- 
tice gases, have been proposed in two- [l9l - l23| and three- 
dimensions [23 - [28j to investigate the thermodynamic 
anomalies presented by water and tetrahedral liquids. 
Some of them were also found to present dynamic anoma- 
lies similar to the ones in liquid water [29l - l3ll |. In these 
models, the distinction between hydrogen bonds and van 
der Waals interactions is essential: these two main in- 
gredients contribute to the appearance of a competi- 
tion between distinct molecular states presenting high 
density (lowly bonded) and low density (highly bonded) 
structures. Nevertheless, in many models, more specific 
molecular interactions, and structures, are used to bias 
the system towards a low density liquid (LDL), at low 
pressures, or a high density liquid (HDL), at high pres- 
sures. As an example, some models use many-body inter- 
actions to unfavor molecular packing in the neighborhood 
of a hydrogen bond [13, [13, ill, 113 ■ Others actually en- 
ergetically favor the LDL states through a repulsive van 
der Waals interaction [2l|, [251 [28l| . Some implement fluc- 
tuating bonding structures [23|, additional unbounded 
molecular states (to stabilize a disordered anomalous liq- 
uid) [13, [H, [13 and some even use ad hoc variations of 
volume with bond formation [23| . 

Considering the increasin g co mplexity found in models 
for water in the literature [32| . simple three- and two- 
dimensional models of liquid water, including only van 
der Waals and hydrogen bond interactions, have been 
investigated with the aim of finding the minimal require- 
ments for water-like anomalous behavior [U [1^, [sO, HH, 
[33 35]. One of these models, the GBHB model proposed 
by Girardi et al.[28j. is a three-dimensional fluid, de- 
fined on a body centered cubic lattice, with first-neighbor 
van der Waals and hydrogen-bond like interactions. It 
presents a phase diagram with two distinct liquid phases 
(high-density and low-density liquid - HDL and LDL) 
besides a GAS phase. Two coexistence lines (GAS-LDL 
and LDL-HDL) ending at critical points were origi nally 
found with Monte Carlo simulations by Girardi et al. [2^ . 
Also, the isobars present temperatures of maximum den- 
sity (TMD) on a line in the pressure-temperature plane, 
resembling qualitatively the scenario emerging from the 
simulations by Poole et al Nevertheless, a qualita- 
tively different phase diagram was found for the same 
model in a recent work by Buzzano and collaborators 
[36| , in which the phase diagrams of a three-dimensional 
model of network forming fiuid ^] were investigated us- 



ing the cluster variational method [37[. With this ap- 
proach they were able to show that the topology of the 
phase diagram of the model was much more complex than 
originally found with Monte Carlo simulations but, at the 
same time, very diverse from the one expected for water. 
It was found that the so-called critical points were indeed 
tricritical points connected to a line of critical points. 
Besides that, another line of critical points was found 



separating the GAS and HDL phases, terminating in a 
critical end point on the GAS-LDL coexistence curve. 

In a more recent paper (38| from the same group, the 
previous analysis was extended by including another two 
three-dimensional models of 'liquid water', also defined 
on the bcc lattice, originalfj^proposed by Bell [13] and 
by Besseling and Lyklema [25]. They revisited the three 
models using the same methodology and the same conclu- 
sion holds for them: in all cases the phase diagrams were 
indeed much more complex than originally expected. In 
the previous analytical studies 24, 25;], phase diagrams 
were oversimplified due to a 'homogeneity' assumption 
on the lattice sites, and by allowing sublattice ordering, 
more stable ordered phases appear and the disordered, 
homogeneous and water-like fluid becomes either unsta- 
ble or metastable (38| . 

Here we investigate a simplified version of a two di- 
mensional associating lattice gas model on the core of 
the Husimi cactus [2l[, considering these recent results 
on lattice models with water-like behavior. The origi- 
nal model was proposed by Henriques and Barbosa and 
studied through Monte Carlo simulations in a series of 
papers [llj, [SO, IH, [13, HI]. In the Henriques-Barbosa 
model each site of a triangular lattice can be occupied by 
a water molecule or empty. A molecule has four bond- 
ing arms (two donors and two acceptors) and two inert 
arms separated by an angle of 180°. All arms lie on 
lattice edges. A HDL was found at low temperatures 
and high pressures for repulsive van der Waals interac- 
tions, while a LDL was found at low temperatures and 
lower pressures. The first Monte Carlo simulations pro- 
vided indications of a coexistence between the HDL and 
the LDL, with the presence of a second critical point 
(SCP) at the end of the HDL-LDL coexistence locus 
[2ll 1331 ]. A temperature of maximum density was also 
found in the neighborhood of this SCP. Variations of this 
model were also investigated through Monte Carlo sim- 
ulations: the distinction between donors and acceptors 
was excluded from the model and distortions were intro- 
duced in the bonding arms [llj. In all cases, the SCP 
and a line of TMD were found to be present, in an in- 
dication of the apparent robustness of these features in 
the phase diagram. More recently, the phase diagram of 
the Henriques-Barbosa model was revisited using simu- 
lations and it was found to be much more complex and 
richer than originally observed [s^ . The new simulations 
suggest that the GAS-LDL coexistence curve ends at a 
tricritical point, and that the LDL-HDL coexistence ends 
at a bicritical point, where the two continuous transition 
lines (GAS-LDL and GAS-HDL) also meet [1|. 

In this work we consider the version of the Henriques- 
Barbosa model without distinction between proton 
donors and acceptors [34]. This simplifying assumption 
does not lead to essential differences in the phase di- 
agrams of this model, particularly with respect to the 
presence of the density anomaly and the HDL-LDL first 
order phase transition The Husimi cactus is built 

with hexagonal plaquettes with a central site (composed 
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by six elementary triangles) as base cells, hereafter called 
hexagons only. This may be seen as a second-order ap- 
proximation on the triangular lattice pol | . Hexagons were 
chosen as a base cells because they are the simplest al- 
ternative we found to reproduce exactly the ground state 
of both ordered phases (LDL and HDL) on the trian- 
gular lattice. We advance that the phase diagram of 
the Henriques-Barbosa model we obtained turned out to 
be very different from the one originally obtained with 
Monte Carlo simulations [s^]- Nevertheless, it is closer 
to the more recent simulations of the model with dis- 
tinction between donor and acceptor arms (soj . In our 
study, the GAS-LDL and LDL-HDL coexistence lines de- 
veloped into two first order phase transitions ending at a 
triple point. In addition to this, a novel first order tran- 
sition line between the GAS and HDL phases appeared 
separating both phases for all pressures. Although our re- 
sults show that the Henriques-Barbosa model may have a 
complex and intriguing phase diagram, in the current for- 
mulation the model seems to be inappropriate for liquid 
water. Nevertheless it does present some water-like fea- 
tures such as a temperature of maximum density (TMD) 
in the fluid phase, which can be used as a starting point 
for more complex two-dimensional models of liquid wa- 
ter. 

In our opinion, in models for complex fluids the com- 
bination of approximate calculations with extensive nu- 
merical simulations are complementary in the study of 
their thermodynamic behavior. Although approximate 
analytical results may be at variance with the correct 
ones for the corresponding model, they may also suggest 
more detailed numerical studies of the model to ascertain 
that the real behavior is found. Besides, it is remarkable 
that a very good agreement was found between the sim- 
ulational and the cluster-variational results for the 3D 
associating lattice gas in [s!]. As will be shown later, 
this is also true for the 2D model studied here. 

This paper is organized as follows. In section [H] the 
model is introduced in more detail on the triangular 
lattice and its ground state is analyzed. We then pro- 
ceed defining the model in a Husimi lattice built with 
hexagons, such that the ground state properties on the 
triangular lattice are exactly reproduced on the Husimi 
lattice. We also present the solution of the model in terms 
of recursion relations and the calculations of the grand- 
canonical potential in the bulk of the tree. In section IIIII 
the thermodynamic properties of the model are studied 
and compared with Monte Carlo simulation data found 
in the literature for the same model. Final discussions 
and the conclusions may be found in section HVl 



II. DEFINITION OF THE MODEL AND 
SOLUTION ON THE HUSIMI LATTICE 

We consider the simplified version of the Henriques- 
Barbosa model on the triangular lattice. Each site of 
the lattice may be either empty or occupied by a single 



molecule. A molecule has four bonding arms, without 
distinction between donors or acceptors of protons, and 
two neutral (non-bonding) arms. The neutral arms form 
an angle 180°, and therefore each particle has three pos- 
sible orientations of the bonding arms. Thus, we are 
considering the symmetric undistorted case discussed in 
[3^ . The possible configurations of a site i will be rep- 
resented by a variable rji, which vanishes if the site is 
empty and assumes the values 1, 2, or 3 if the site is oc- 
cupied in one of the possible orientations of the bonding 
arms. Repulsive van der Waals interactions e > exist 
between particles on first neighbor sites, and an energy 
7 < corresponds to each hydrogen bond on the lattice. 
Therefore, if I7I > e, a pair of particles on first neigh- 
bor sites with an hydrogen bond between them is asso- 
ciated to a net negative energy and thus the interaction 
becomes attractive. Since we will study the model in the 
grand-canonical ensemble, an activity z — e:xjp{fj,/kBT) 
corresponds to each particle on the lattice, where ^ is 
the chemical potential. We may relate the parameters 
used here and those chosen in reference [s^l , there a pair 
of first-neighbor sites occupied by particles with a hy- 
drogen bond between them corresponds to an energy —v 
and if no hydrogen bond is present this energy is —v + 2u 
(4l| . Therefore, we have e = —v + 2u and 7 = —2u, and 
we notice that for u/v — I we have \ — 2. Since the 
simulations in f34| were done for this particular choice, 
we restrict our numerical calculations to this particular 
case. 

Three phases were found in the ground state in ear- 
lier investigations [21], [sJl : The GAS phase corresponds 
to the empty lattice, and is stable at low values of the 
chemical potential; as the chemical potential is increased, 
the low-density liquid (LDL) becomes stable, in which a 
fraction p = 3/4 of the sites are occupied by particles 
and all lattice edges between two particles are occupied 
by hydrogen bonds. For still higher chemical potentials, 
a high-density liquid (HDL) becomes stable, in which all 
sites are occupied and therefore p — I- In Fig. [T]both 
liquid phases in the ground state are depicted. 

To describe the sublattice structure of the LDL phase 
it is necessary to introduce four sublattices: three of 
them formed by sites on the borders of the hexagons 
and the central sites (see Fig. [5] (a)). This immedi- 
ately leads to a fourfold degeneracy (corresponding to 
the placement of a hole in four different sites), which is 
captured by the Husimi tree if the homogeneity assump- 
tion is avoided ([1^). The sublattices are equivalent in 
the HDL but we take care of the orientational order. This 
leads to a threefold degeneracy (corresponding to possi- 
ble orientations for water molecules in the lattice), and 
this behavior is also captured by the tree without the 
homogeneity assumption. As will be shown later, the 
states of the central site of each hexagonal plaquette will 
always be summed in the recursive equations obtained 
from the hierarchical structure of the lattice. Consider- 
ing this (and also to simplify our notation), we use only 
the three sublattices for the sites on the perimeters of the 
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FIG. 1: A representation of the two ordered phases of the 
model with the configuration of each particle identified by 
the orientation of the inert arms (upper panel). In both high 
density liquid (middle panel) and low density liquid phases 
(bottom panel) hydrogen bonds are indicated using full lines 
lattice edges with van der Waals interactions only are drawn 
with dashed lines, while the other lattice edges are represented 
by dotted lines. 



hexagons, as is shown in Fig. [S] 

Let us now discuss the ground state of the model in 
some detail. In the GAS phase all sites are empty and 
we will associate a vanishing energy to this configura- 
tion, Egas = 0. The LDL phase on the Husimi lattice 
is fourfold degenerate, characterized by empty sites ei- 
ther at the center of each hexagon or at the sites of one 
of the three sublattices A, B, and C. Recalling that all 
edges between first neighbor sites occupied by particles 
have hydrogen bonds on them, the energy per hexagon 
(including the chemical potential term) will be: 

ELDL=Q{e + l)-i^J^. (1) 

where we remember that each particle on the vertices of 
the hexagons is shared by two plaquettes. In the HDL 
phase, all sites are occupied and 8 of the 12 edges of each 
hexagon are occupied by hydrogen bonds, while the re- 
maining 4 are not. Thus, there are three possible config- 
urations of the hydrogen bonds. The energy per hexagon 
in this phase is; 

Ehdl = 12e + 87 - 4^. (2) 



It is easy to find which phase corresponds to the min- 
imum energy for given parameters e, 7, and ^. Using 
the vdW interaction e as the energy scale, we may de- 
fine the dimensionless variables 7 = and p, = /x/e. 
The ground state corresponds to the GAS phase if /2 < 
2(1-7), to the LDL phase if 2(1-7) < M < 2(8-7), and 
to the HDL if /2 > 2(3 — 7). As observed above, these val- 
ues are the same as the ones found for the ground state 
on the triangular lattice [33 |. 

A. Recurrence relations on the Husimi cactus 

As usual, we start defining partial partition functions 
for rooted subtrees, fixing the configuration of the root. 
One of these subtrees is shown in Fig. [2] There are 3x4 
configurations of the root sites, so we define 12 partial 
partition functions gi, i = 1,2,. ..,12, where the index i 
stands for the root site configuration. We may associate 
the configurations (s,?]), where s = A,B,C stands for 
the sublattice and 77 = 0, 1, 2, 3 for the site configuration, 
to the indices i of the partial partition functions following 
the convention indicated in table ID It is also useful to 
define the configurations of the bonding arms of a par- 
ticle in a way which may be applied to all sites of the 
tree (only sites at the perimeter of the hexagons are con- 
sidered, since we will sum over the configurations of the 
central sites). We therefore consider a particular site with 
a particle on it, which will belong to two hexagons in dif- 
ferent generations of the tree, and imagine that we circle 
around the site clockwise, starting outside the hexagons. 
The configuration variable 77 associated to this site will 
be equal to the number of lattice edges we cross until 
the one where one of the inert arms of the particle are 
located is reached, added with one. This definition is 
illustrated in Fig. [21 We proceed considering the opera- 
tion of attaching 5 subtrees with M generations to a new 
root hexagon, building a subtree with M + \ generations. 
Summing over the 4^ = 1024 possible configurations of 
the root sites of the M-generations subtree, we will arrive 
to recursion relations for the partial partition functions, 
which are of the form: 

1024 / 3 \ 12 

5f = E E^"^"^p"^"-'"-^" n(ff )"'^■■^ (3) 

where rij^k, Pi,j,k, and are the number of parti- 

cles, number of pairs of particles in first neighbor sites 
and number of hydrogen bonds for each contribution j 
to the partial partition function of gf^~^^, given that 
the configuration 77 of the central site is equal to k. 
ujp — exp[— e/(fcBT)] and Wf, — exp[— 7/(fcBr)] are the 
Boltzmann factors associated to the van der Waals inter- 
actions and hydrogen bonds, respectively. We notice that 
the activity of the particle which eventually is placed on 
the root site is not considered at this level. The expo- 
nents ei,j^i assume integer values between and 2. We 
remark that for each pair of indices (i, j) there will be at 
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FIG. 2: a) Definition of the sublattices. b) A subtree with three generations, c) Definition of ri for a site occupied by a particle. 
The numbers on the lattice edges correspond to the values of the variable 77 if the inert bonds of the particle are placed on 
these edges. 





i 


i 


{A,l), (A 2), (A, 3) 


1, 2, 3 


10 


(3,1), (5,2), (5,3) 


4, 5, 6 


11 


(C,l), (C,2), (C,3) 


7, 8, 9 


12 


{A,0), (5,0), (CO) 


10, 11, 12 





TABLE I: To each possible state of a site, specified by its 
sublattice s and configuration rj (first column), an index i is 
assigned (second column). The indexes i which appear in the 
denominator of eq. ((4]) are shown on the third column. 



most five nonzero exponents e,;.j^^, since this is the num- 
ber of subtrees with M generations linked to the new 
hexagon at the root. The exponents Sij^i depend of the 
index i only because the number of incident subtrees with 
root sites in each sublattices depends of the sublattice of 
the root site of the new subtree. 

In similar calculations, often the recursion relations 
are obtained by hand, usually summing the contributions 
with some graphical aid. In the present case, due to the 
large number of contributions, this procedure is very te- 
dious and therefore errors are quite frequent. Although 
we actually obtained the recursion relations explicitly in 
this way, using symmetries to generate the expressions 
for the recursion relations, these expressions are much 



too large to be given here. To assure that the recursion 
relations are free of errors, we also decided to write a 
rather simple code which generates the sets of 24 integer 
numbers Uj^k, Pi,j,k, bij^k and Cij^i for each contribu- 
tion j to the recursion relation for g^^^^ , similar to what 
was done by Zara and Pretti in a model for RNA on the 
Husimi lattice [42| . Since we are interested in the behav- 
ior of the model in the thermodynamic limit, we should 
consider fixed points of these recursion relations. As ex- 
pected, however, the partial partition functions diverge 
in this limit. So, we may define ratios of these functions, 
which may approach a finite value as M — ^ 00. We thus 
define the ratios Ri, i = 1,2,...,9, dividing each partial 
partition function with a particle at the root site by the 
partial partition function with an empty root site in the 
same sublattice. This leads us to the ratios Ri — gi/gj, 
where the values of i are shown in table HI We may then 
obtain recursion relations for the ratios from the ones for 
the partial partition functions, Eq. |31 They are: 



(4) 
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B. Densities in the core of the Husimi cactus 

In order to obtain densities in the central region of the 
tree, we consider the operation of attaching 6 subtrees to 
the central hexagon, which leads to an expression for the 
partition function of the whole tree: 
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(5) 



e=i 



where Nj^k, Pj,k and Bj^k are the total number of parti- 
cles, nearest neighbors and hydrogen bonds on the central 
hexagon with the configuration of the border sites given 
by j and the central site in the configuration k, and Ejj 
is the number of border sites with configuration £. 

Again the set of integer exponents was generated by a 
computer program, as well as manually, and both proce- 
dures lead to the same final results. Now, for example, 
the density of particles (defined here as the number of 
particles divided by the number of sites) in the central 
hexagon will be given by: 

z dYM 

where p is in the range [0, 1]. We notice that the activities 
of the sites at the perimeter and the center of the cen- 
tral hexagon of the tree are considered in expression ([5]) , 
so the factor 7 assures the proper normalization of the 
density. In other words, the numbers Nj^k in expression 
^ are in the range [0,7]. A similar procedure leads to 
expressions for the densities of hydrogen bonds and van 
der Waals interactions per site. Dividing both the nu- 
merator and the denominator of the expressions for the 
densities by {gfl gff 5*1)^ ^lay express them in terms 
of the ratios and the parameters of the model. Thus, for 
example: 



1 2^j=i yz^k=o 



M\Ei 



7 Y^4096 fsr^3 



(7) 

To obtain the thermodynamic behavior of the model, 
we may iterate the recursion relations until a fixed point 
for the ratios is reached with the required numerical 
precision, and then calculate the densities at the center 
of the tree. The convergence of the recursion relations 
generally is quite fast. In certain regions of the parameter 
space, more than one fixed point may be stable, signaling 
coexistence of phases. To locate the first order transition 
in such cases it is necessary to compare free energies of 
different phases. An expression for the grand-canonical 
free energy is obtained in what follows. 



C. Grand-canonical free energy 

To obtain the grand-canonical free energy of the model 
in the core of the tree, we may proceed following the pre- 
scription proposed by Gujrati [40]. For this purpose, it 



is convenient to notice that if we connect the central site 
of each hexagon to the central sites of the first-neighbor 
hexagons, we end up with a Cayley tree with coordina- 
tion q — 6 and ramification a = q — 1 = 5. Now we 
may assume that the total free energy of the tree is the 
sum of the free energies associated to each hexagon. This 
takes care of the sublattice structure of the model. The 
hexagons may then be classified in generations identified 
by the index m, starting with the ones placed on the 
surface, for which m — 1 and ending at the central site 
(m = M for a tree with M generations of hexagons). It is 
natural, considering the structure of the tree, to assume 
a radial symmetry for the local free energies, so that we 
will represent by 0m (w) the free energy of a hexagon in 
the m'th generation of a tree with a total of M genera- 
tions. The total free energy of the tree will then be: 



M 



= ^ NM(m)<j>M{ni), 



m—1 



where Nm{M) = 1 and 

NM{m) = q<j^^-"'-\ m = 1, 2, . . . , M - 1 



(8) 



(9) 



are the numbers of hexagons in the m'th generation. Now 
we may see that: 

$Af+l - G^M ^(i)M+l{M + l)-G(j)M{M) + q<j)M+l{M) 



M 



-9 5]a^^-™[(/.M+i(m)-(/.M(m)]. 



(10) 



It is now reasonable to assume that in the thermody- 
namic limit we should have the free energies per hexagon 
in the core of the tree approaching a limiting bulk value 
so that (j)M+i{M + 1) = 0m(M) = (j)M+iiM) = (j)b 
for M ^ oo. Close to the surface, the free energies per 
hexagon should be functions of m, but we may assume 
4>M+i{iTi) — 4>M{m) — > when M oo. Therefore, we 
may conclude from equation (jlOp that: 



1 



($ 



M) 



(11) 



in the thermodynamic limit Af — > oo and for the model 
we are considering here cr = 5. Actually, this expression 
may also be obtained from a somewhat stronger assump- 
tion that the free energies per hexagon in the thermo- 
dynamic limit should assume only two values: (pb on the 
surface m — 1 and 0f, in all other cases Also, we no- 
tice that the original argument for the calculation of the 
free energy was also recently generalized for Husimi trees 
with a sublattice structure in Semerianov and Gujrati 
[31 ■ The derivation presented here is different from the 
ones originally proposed in [i^l and Q, but the results 
are the same |45l |. 

Since the partition function on a M-generations tree is 
Ym, we have: 



1 



-ksTln 



Ym+i 

^ M 



(12) 
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Substituting the partition function ([5]) in this expression, and expressing the sums in terms of the ratios, we obtain: 



Ym+1 



(13) 



Now we may use the recursion relations Eqs. (jSj to express the partial partition functions for subtrees with M + 1 
generations in terms of the ones with M generations, and finally will arrive at the expression for the second fraction 
in expression (jl3p 



{gT'gli^'gr.^'f 
{g'llgfigfiY' 
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n 




71 A I- Pi,i,k ^i, 7, A: 



r 



(14) 



Therefore, we see that we may express the bulk free 
energy per hexagon as a function of the parameters of 
the model and the ratios and in the thermodynamical 
limit it will converge to a fixed point value. Finally, since 
we are in the grand-canonical ensemble, we have that the 
pressure is P — — where $ is the grand-canonical 
potential and V the volume. Associating a volume vq to 
each site of the lattice and recognizing (^t as the grand- 
canonical potential per hexagon for the solution on the 
Husimi tree, the pressure may be written as 



P = -0fc/4i;o, 



(15) 



where it should be stressed that we have four sites per 
hexagon in the core of the tree. 



III. THERMODYNAMICAL BEHAVIOR OF 
THE MODEL 

To study the thermodynamical behavior of the model 
on the Husimi lattice, we define reduced intensive or 
fieldlike thermodynamic variables (temperature, pressure 
and chemical potential) as T = /cgT/e, /I — /i/e and 
P = Pva/e. Considering expressions (fT2l) and (fT5|) . the 
reduced pressure is given by: 



In 



P^T- 



8 



(16) 



As mentioned before, we limited our study on the par- 
ticular case 7 = 2, for which MC simulations were found 
in the literature. For fixed values of T and /l, we iter- 
ate the recursion relations ^ for the ratios of partial 
partition functions, and once the fixed point is reached, 
we determine the mean numbers of particles, hydrogen 
bonds and van der Waals interactions per lattice site, 
which are represented by p, Vhb, and v^y/i respectively. 
The densities of hydrogen bonds and van der Waals in- 
teractions per site, normalized to be in the range [0, 1], 



will be phb = vhb/'^ and p^w = v^w /"i- Finally, the 
pressure may be also obtained at the fixed point. 

The phase diagram of the model was found using this 
procedure, being presented in the (T, P) plane on Fig. 
[3l The three phases used in our ground state analysis 
were also found at finite temperature and coexistence 
lines between these phases were calculated by requiring 
the identity of their bulk free energies. All transitions are 
discontinuous and a triple point, located at P = 2.997, 
T = 0.835, and Jl = 1.959, was found. The coexistence 
lines meeting at the triple point satisfy the thermody- 
namical requirements for this situation, such as the 180 
degree rule [46j . 

Our phase diagrams are qualitatively different from the 
one originally suggested based on Monte Carlo simula- 
tions [34|, where the two coexistence lines start at low 
temperatures and end at critical points. They are also 
different from the cluster variational results obtained for 
several waterlike models on the hcc lattice[3^ H^, where 
the coexistence lines end at tricritical points. In Fig.[3](c) 
the MC simulation results presented in [sl] are also 
shown and a good agreement on the location of the phase 
transitions is found between those and our results, except 
on the low-temperature region of the LDL-HDL coexis- 
tence line. Along this region the simulations present a 
rather large positive slope, and since the ground-state 
value is exactly known (P = 3), a large negative slope 
should occur in the LDL-HDL coexistence line at low 
temperatures. This effect is even enhanced if we recall 
that, due to the third law of thermodynamics and the 
Clausius-Clapeyron equation, the coexistence curve has 
to be horizontal at vanishing temperature (47j . A simi- 
lar situation was found in the simulations of the model 
with distinction between donor and acceptor arms [2]] |. 
where this point is discussed, particularly with respect 
to the implications of the Clausius-Clapeyron relation. 
Although these simulations have been recently revisited 
[39} , the new results for the LDL-HDL coexistence line do 
not include temperatures low enough to reach the region 
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FIG. 3: (color on line) a) Phase diagram {f x P) of the 
Henriques-Barbosa model on the Husimi lattice, as defined on 
section [IT] Dashed lines are discontinuous transitions which 
meet at triple point represented by a full circle (red on-line). 
The full lines are the stability limits of the GAS (no symbol, 
red on line), LDL (squares, blue on line), and HDL (trian- 
gles, green on line) phases. The dotted and dash-dotted lines 
(black on line) are the TmD and TMD, respectively, b) A de- 
tail of Fig. a) where the lines of density anomalies are more 
visible, c) Present results for the coexistence lines and triple 
point are compared to the first order phase transitions (cir- 
cles) and TMD (triangles) from Monte Carlo simulations of 
Balladares et al [341 ]. 



we are discussing here. In our calculation, the HDL-LDL 
coexistence curve starts with zero slope at vanishing tem- 
perature. As the temperature increases, the slope has a 
small positive value, then the curve presents a maximum 
and the slope becomes negative close to the triple point. 
These features are not visible in the scale of Fig. |31 It is 
interesting to notice that the estimated location for the 
LDL-HDL critical point in the simulations is quite close 
to the triple point in our solution. 

We carefully verified if the transitions are actually dis- 
continuous by studying the stability limits of the fixed 
points associated to each phase. These limits may be 
found calculating the jacobian of the recursion rela- 
tions (gl) 



dR 



M+l 



(17) 



at the fixed point (M oo), and then requiring the 
absolute value of the largest eigenvalue of the jacobian 
to be equal to one. In Fig. ^a.) the stability limits of 
all phases are shown. Although in part of the GAS-LDL 
coexistence line the stability limit of the LDL phase is 
very close to the transition, they are never coincident, 
thus assuring the discontinuity of the transition. 

In order to find out if the stability limits of the fixed 
points are in fact the spinodals (thermodynamic stability 
limits), we also calculated the eigenvalues of the hessian 
associated with the phases. For the ordered phases (LDL 
and HDL) we found a good numerical coincidence of the 
spinodals and the stability limits everywhere. For the 
GAS phase, at low temperatures, we were able to assure 
numerically the coincidence between these curves, but 
at higher temperatures we had numerical problems to 
evaluate the elements of the hessian, which are second 
derivatives of the potential. 

An interesting point is that the LDL-GAS coexistence 
line has two regions with slopes of different signs, showing 
a reentrant behavior. The p, x T phase diagram is quite 
similar to the P xT phase diagram shown in Fig. [3] The 
change of the sign happens at a point which is located at 
fimax = 0.137, Tmax = 1.029, and Pmax = 1.686. Since 
the GAS phase has a larger entropy at the coexistence 
with the LDL phase, the Clausius-Clapeyron relation in- 
dicates that the particle density should be lower for the 
GAS phase than for the LDL phase in the part of the 
coexistence curve with pressures lower than Pmax- At 
{Tmax, Pmax) ^ the densities of both phases are identical 
and in the remainder of the coexistence the density of 
the GAS phase is higher. In fact, the equal densities at 
this point are confirmed in Fig. SJa), where the tem- 
perature is shown as a function of the density of parti- 
cles at coexistence. It is important to remind that GAS 
and LDL phases are not identical on this point (in this 
case it would be a critical point). As can be observed in 
the phase diagram with the density of hydrogen bonds, 
instead of the particle density, shown in Fig. Sfb). Al- 
though not presented here, densities of vdW interactions 
are also different for both phases on this point. 
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FIG. 4: (color on line) Temperature x density diagrams, a) 
Particle density, b) Hydrogen bond density. The temperature 
of the triple point is indicated by a dashed line. Densities of 
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line) corresponds to the HDL phase, the dashed line (blue on 
line) to the GAS phase and the dot-dashed line (red on hue) 
the the LDL phase. The dotted horizontal line corresponds 
to the three-phase coexistence (triple point). 



Another relevant question is the location of the points 
of maximum density in the isobars. We found that iso- 
bars for the densities of particles as functions of the tem- 
perature do actually present a maximum at pressures 
above Pmax- This TMD is located in the metastable 
extension of the GAS phase inside the LDL phase, in 
Figs.|3Ja) andinjb), we represented the location of these 
metastable TMD of the GAS with the dash dotted line. 
It ends at the point of maximum temperature of the GAS 
spinodal, as may be seen in the detail (Fig. [SUb)). We 
also found a temperature of minimum density line (TmD) 
inside the LDL phase, shown in Figs, ^a) and (b) as a 
dotted line. This line occurs also at pressures higher 
than Pmax and, unlike the TMD, which is located in the 
metastable GAS phase, the TmD line covers both stable 
and metastable regions of the LDL, ending at the point 
of maximum temperature of the LDL spinodal, a detail 



also more visible in Fig. |3l^b). It is actually expected that 
lines of vanishing thermal expansion coefficient should 
end at the points where the corresponding spinodals 
change the sign of their slope [1^. Finally, it should be 
mentioned that similar findings were reported by Pretti 
and Buzzano in their homogeneous cluster variational 
study of the symmetric Roberts-Debenedetti model ^3] ■ 

In Fig. [5]we show some isobars for the densities of par- 
ticles and hydrogen bonds. The results of the present 
calculations are represented by the broken lines, and the 
symbols close to the isobars are the MC simulation results 
obtained by Balladares et al [13] . We notice a good quan- 
titative agreement between them and our calculations, at 
least not too close to the coexistence line. In the P x T 
diagram shown in Fig. [3jc), the locations of the TMD 
points found in the simulations presented in (33 | are rep- 
resented by triangles, and in general we may notice that 
they are located at temperatures larger than the ones of 
coexistence. These estimates actually correspond to the 
maxima in the density at the coexistence curve, and the 
fact that they lie above the coexistence curve may be due 
to finite-size effects. As another possibility, the Bethe 
lattice approximation introduced here could underesti- 
mate the location of the TMD due to the absence of a 
LDL-GAS critical transition, observed in simulations, at 
least for the model with distinction between donors and 
acceptors [39j . In principle, the presence of such a criti- 
cal line could increase the entropy-volume cross fluctua- 
tions and shift the TMD line {VksTa = {SV5S) = 0) to 
higher temperatures. For example, this kind of TMD un- 
derestimation happened in Bethe lattice solution of the 
Bell-Lavis model of liquid water fssj , when compared to 
Monte Carlo simulations [49]. 

In Fig. in^b) the mean number of hydrogen bonds per 
particle {nHB = vhb/p = '^Phb/p) is depicted at con- 
stant pressure as a function of the temperature. Again, a 
good quantitative agreement between our results and the 
MC simulations was found. In the simulations, crossing 
of different isobars was reported and its physical origin 
was discussed [S^l, in relation to density anomaly. Nev- 
ertheless, our calculations suggest the possibility of the 
crossing being a consequence of the discontinuous tran- 
sition between the LDL and GAS phases and finite-size 
rounding effects in the isobars. 



IV. FINAL DISCUSSION AND CONCLUSION 

In this paper we solved the Henriques-Barbosa model 
with symmetric arms [34] on a Husimi lattice built with 
hexagons. Two liquid and one gas phase are present in 
the phase diagram, but qualitative differences are found 
when compared with the phase diagram which was ob- 
tained with simulations. All transitions we found are 
discontinuous, and also a triple point was found where 
all phases coexist with different densities. We carefully 
checked if the transitions are really discontinuous, since 
in part of the coexistence loci the discontinuities are 
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rather small. Thus we assured that the stability lim- 
its of the fixed points associated to the coexisting phases 
are never coincident. Also, it may be seen in Figs. |4]that 
the densities present a discontinuity at the coexistence 
line, although it may be rather small, particularly in the 
neighborhood of the point of maximum temperature in 
the GAS-LDL coexistence line. 

Recently, more detailed simulations were reported on 
this model with distinction between donor and acceptor 
arms ^39] , and a diagram closer to the one we present here 
was found. The difi^erence is that the GAS-LDL transi- 
tion line is discontinuous at low temperatures, but be- 
comes a critical line when the temperature is increased, 
thus a tricritical point is present. Also, the HDL-GAS 
transition appears to be continuous in the new simula- 
tions. The LDL-HDL line is always discontinuous, so 
that in the simulations the point which corresponds to 
the triple point in our phase diagrams appears as a bicrit- 
ical point (which was called wrongly as a tricritical point 
in the caption of Fig. 3 in reference [s^). Nevertheless, 
a direct comparison between the present calculation and 
these new simulations may not be done, since the dis- 
tinction between donor and acceptor arms leads to an 
increase of the entropy of the model. It is not impossible 
that a transition found to be discontinuous in mean-field 
like approximations turns out to be continuous in simu- 
lations or more precise calculations, such as series expan- 
sions. However, we notice that the results of the calcu- 
lations presented here show, in general, good agreement 
with data furnished by simulations, as was also noticed 
by Buzano and collaborators in their study of an asso- 
ciating lattice gas model, with tetragonal symmetry, on 
the bcc lattice [36] . It is worth mentioning that the phase 
behavior presented by the most recent simulations of the 
Henriques-Barbosa model is different from that found for 
the associating lattice gas model studied with the clus- 
ter variational method on Ref. [S^. There, instead of a 
bicritical point, a tricritical and a critical endpoint are 



present. We are presently studying this model with the 
same methods implemented on this paper. 

Although the results presented here for the coexistence 
lines in the pressure-temperature phase diagram agreed 
well with the data of simulations in jsj], this is not the 
case for the low-temperature region of the LDL-HDL co- 
existence curve, where the simulations suggest a mini- 
mum while a smooth behavior was found in our results. 
The interpretation of this apparent minimum in relation 
with the Clausius-Clapeyron equation seems unclear to 
us, and we believe on the possibility that these results 
might be spurious. Possibly, longer equilibration times 
for MC simulations should be considered on this low- 
temperature region. 

If we adopt the qualitative phase diagram which 
emerges from our calculations, the TMD found in the 
simulations would correspond to the coexistence line. 
Nevertheless, it is also possible that the absence of a crit- 
ical line results in an overall decrease of the temperatures 
of the TMD line, as in the case of the Bell-Lavis model 
(35I |49| . The crossing of isobaric curves for the density 
of hydrogen bonds as function of the temperature was 
observed in the simulations, but here it may be seen as 
a consequence of rounding finite size effects for the dis- 
continuous transition at the LDL-GAS coexistence curve, 
with no relation to the TMD. The LDL-GAS coexistence 
curve actually is a very weak discontinuous transition in 
the region where the simulations suggested the presence 
of a critical point. This indicates that the question of the 
order of the transition in this region should be studied 
very carefully in simulations. 

The results presented here support that the phase di- 
agram of the symmetric Henriques-Barbosa model does 
not present a second critical point, as originally suggested 
through Monte Carlo simulations [s^l- Our theoretical 
phase diagram is actually closer to the results of more 
recent simulation of a original (asymetrical) model {s^: 
phase transitions are found in both cases, but second or- 
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der transitions are not present in our theoretical results, 
this may be due to the hmitation of the Husimi lattice 
solution to capture long-range correlations. 

Considering the technical aspects of the Husimi lat- 
tice, our calculation shows the importance of a careful 
choice of the sublattice structure to capture the correct 
phase diagram of not so simple lattice models. A sim- 
pler homogeneous lattice structure would lead the sys- 
tem to thermodynamic states which may be unstable or 
metastable in the more general parameter space used here 
(results not shown). On this sense, it is worth mention- 
ing that usually the simplest choice for the plaquette in 
Husimi lattice calculations to approximate the behavior 
of models on regular lattices is the elementary polygon in 
the original lattice. Therefore, the simplest choice in the 
present case would be a Husimi lattice with triangular 
plaquettes and coordination number equal to 6. How- 
ever, such a choice would probably not lead to results 
comparable to the ones obtained in simulations. The 
rather unusual choice of plaquettes we adopted here is 
certainly responsible for the good quantitative agreement 
between our results and the simulations. This technical 
discussion is extremely relevant to the subject since one 
can find in the literature lattice models presenting water 
like behavior (including thermodynamic anomalies and 
interesting phase diagrams) which were studied without 
a more detailed sublattice analysis [s^, HH . 

We finish remarking that many lattice models pro- 
posed to investigate waterlike anomalous behavior were 
found to present phase diagrams which were more com- 
plex than originally expected. When sublattice ordering 



is properly considered 



even the simplest models do 



present at least a single critical line, and many of them do 
not present gas-liquid phase transition ending in a critical 
point or a temperature of maximum density in a stable 
disordered fluid without a sublattice structure. The ar- 
bitrary use of the homogeneneity assumption can result 
in some interesting phase diagrams but, in our opinion, 
this assumption is an artifact which eventually may hide 
the instability of the homogeneous phase at low temper- 
atures. Thus, this homogeneous solution could not be 
considered as a true implementation for a random lattice 
useful for representing fluids, as recently proposed [ssj . 
Considering the results obtained here and in other re- 
cent papers [H, HI] , the issue of finding a simple lattice 
model with minimal waterlike behavior seems to be far 
from being resolved. Monte Carlo simulations are cer- 
tainly needed to determinate the 'exact' phase diagram 
but it seems that a modeling breakthrough will be needed 
to achieve a better description of liquid water. 
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